Pressure-induced structural and dynamical changes in liquid Si. An ab initio study 
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The static and dynamic properties of liquid Si at high- pressure have been studied using the 
orbital free ab-initio molecular dynamics method. Four thermodynamic states at pressures of 4, 
8, 14 and 24 GPa are considered, for which X-ray scattering data are available. The calculated 
static structure shows qualitative agreement with the available experimental data. We analize the 
remarkable structural changes occurring between 8 and 14 GPa along with its reflection into several 
dynamic properties. 
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I. INTRODUCTION. 

The intriguing properties of Silicon along 
with its technological importance have stim- 
ulated intensive theoretica l 1 ' 2 ' 3 ' 4 ' 5 ' 6 '?'?' 9 and 
experimental 10 ' 11 ' 12 ' 13 ' 14 ' 15 ' 16 work. Its high-density 
forms include the semiconducting and covalent crys- 
talline and amorphous phases and the metallic liquid 
phase. Upon melting it undergoes a semiconductor-metal 
transition, a density increase of f=a 10% and significant 
changes in the local structure which evolves from an 
open one, with a tetrahedral fourfold coordination, to 
a liquid structure with w sixfold coordination. In crys- 
talline Si (c-Si) the semiconducting diamond structure 
contracts with pressure and transforms at 12 GPaii to 
the metallic white-tin structure and then to the metallic 
simple hexagonal structure at 16 GPai£. The local 
structure of liquid Si (1-Si) at the triple point (TP) is 
somewhat similar to high-pressure forms of c-Si, and it 
has been suggested that 1-Si might consist of a mixture 
of diamond-type and white-tin-type structures with the 
proportion of the latter increasing with pressure. 

Within this backdrop, Funamori and Tsuji 1 ^ have re- 
cently carried out X-ray (XR) diffraction experiments to 
determine the static structure of 1-Si at pressures of 4, 8, 
14 and 23 GPa and temperatures about 50 K above the 
melting point at the pressure. From an analysis of the 
static structure factors S(q) and the associated pair dis- 
tribution functions g(r), Funamori and Tsuji 1 ^ concluded 
that 1-Si up to 8 GPa has a local structure intermediate 
between the diamond-type and the white-tin-type. But, 
between 8 and 14 GPa drastic structural changes were 
noted with 1-Si transforming to a denser structure simi- 
lar to that of 1-Sn at ambient pressure as evinced by the 
strong similarities between the S(q)'s of 1-Sn and 1-Si at 
14 GPa. 

Prompted by these experimental developments, we 
have performed an ab-initio molecular dynamics (AIMD) 
study of several static and dynamic properties of com- 
pressed 1-Si at the thermodynamic states addressed by 
Funamori and Tsuji^S. Of particular interest is the re- 
flection of the reported structural changes in the dynamic 
properties. Our AIMD method is based on density func- 



tional theorjiia (DFT) which, for given nuclear positions, 
allows the calculation of the ground state electronic en- 
ergy and yields the forces on the nuclei via the Hellmann- 
Feynman theorem, enabling MD simulations in which 
the nuclear positions evolve according to classical me- 
chanics whereas the electronic subsystem follows adia- 
batically. Most AIMD methods are based on the Kohn- 
Sham (KS) form of DFT (KS-AIMD methods) which 
treats the electron kinetic energy exactly, but which at 
present, poses heavy computational demands limiting the 
size of the systems to be studied as well as the simulation 
times. Some of these constraints can be relaxed by the 
so-called orbital-free ab-initio molecular dynamics (OF- 
AIMD) method, which approximates the electron kinetic 
energy but disposes of the electronic orbitals of the KS 
formulation. The method allows simulations in which 
the number of variables describing the electronic state is 
greatly reduced so that larger samples (several hundreds 
of particles) can be studied for longer simulation times 
(tens of ps). 

Theoretical studies of 1-Si have mainly focused on 
static structural properties for thermodynamic states 
near the TP. Most studies were classical MD simula- 
tions using effective interatomic potentials constructed 
either empirically by fitting to experimental datai^ or 
derived from some approximate theoretical model^ 4 * 5 -. 
Recently, KS-AIMD calculations^^ have been reported 
which address electronic and static properties. Stich et 
a£ and Chelikowsky et a& have reported KS-AIMD cal- 
culations for 1-Si for 64 particles, using non-local pseu- 
dopotentials and the local density approximation. A sub- 
sequent calculation^ used 350 particles and an improved 
treatment of electron exchange and correlation. Recently, 
we have carried out an OF-AIMD simulation 9 for 1-Si 
near the TP for 2000 particles using a first principles 
local pseudopotential. Both static and dynamic proper- 
ties were calculated with results in good agreement with 
the available experimental data, supporting the validity 
of the OF-AIMD for treating systems such as 1-Si which 
show some remnants of covalent bonding and are not fully 
metallic. 

On the experimental side, besides the aforementioned 
XR experiments of Funamori and Tsuji 1 ^, we also quote 
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the availability of both neutron scattering (NS) 10 and X- 
ray (XR)AAii2iiii diffraction data as well as the recent in- 
elastic X-ray scattering (IXS) data of Hosokawa et oA^i^ 
which have provided information on the dynamic struc- 
ture of 1-Si near TP. 

In the next section the orbital-free ab-initio molec- 
ular dynamics (OF-AIMD) scheme is described briefly 
with emphasis on the electronic kinetic energy functional 
and the local pseudopotential used to characterize the 
electron-ion interaction. In section II I II the results of the 
ab initio simulations for several static and dynamic prop- 
erties are presented and compared with the available 
experimental data. Finally, conclusions are drawn and 
ideas for further improvements are suggested. 



II. THEORY. 

The OF-AIMD method used in this study is described 
fully in earlier workSS, and has previously been used to 
study 1-Si near the TP 9 . In summary an explicit density 
functional for the electronic energy is minimized itera- 
tively for each ion configuration, the forces on the ions 
are found using the Hellman-Feynman theorem and the 
ion positions and velocities are updated by solving New- 
ton's equations. The approximate electron kinetic energy 
functional which correctly gives the Thomas-Fermi and 
linear response limits is based on the von Weizsacker term 
plus a correction which uses an averaged densitjiSS. The 
local electron-ion pseudopotential was constructed, for 
each thermodynamic state, according to the procedure 
described in reference^. 

Simulations have been performed for 1-Si in the four 
thermodynamic states listed in Table [I] These corre- 
spond to pressures of 4, 8, 14 and 23 GPa and tem- 
peratures about 50 K above the melting point for each 
pressure^. Each simulation used 2000 ions in a cubic 
cell with periodic boundary conditions and size appro- 
priate for the ionic number density, pi. The square root 
of the electron density was expanded in plane waves up 
to a cutoff energy Ec nt = 15.75 Ryd. The Verlet leapfrog 
algorithm with a timestep of 3.5 x 10~ 3 ps was used to 
update the ion positions and velocities. Equilibration 
lasted 10 ps. and the calculation of properties was made 
averaging over a further 65 ps. For comparison, we men- 
tion that the KS-AIMD simulations for 1-Si near the TP 
lasted 1.2 psA 0.9 psi and 1.0 ps£, which precludes its 
application to the study of most dynamical properties. 

Several liquid static properties were evaluated during 
the simulation: pair distribution function, static struc- 
ture factor and bond angle distribution, as well as vari- 
ous dynamic properties, both single-particle ones: veloc- 
ity autocorrelation function, mean square displacement, 
and collective ones: intermediate scattering functions, 
dynamic structure factors, longitudinal and transverse 
currents. The calculation of the time correlation func- 
tions (CF) was performed by taking time origins every 
five time steps. Several CF are also dependent on the 



wave vector q = | q | . 

TABLE I: Input data for the different thermodynamic states 
studied in this work, pi is the total ionic number density and 
T is the temperature which have been taken from Ref^-. 



P(GPa) 


Pi (A" 3 ) 


T (°K) 


4 


0.058 


1503 


8 


0.060 


1253 


14 


0.067 


1093 


23 


0.071 


1270 



III. RESULTS 
A. Static properties. 

The simulations yield directly the pair distribution 
function, g(r), and the static structure factor S(q). Fig- 
ure Q] shows the calculated S(q)'s along with the corre- 
sponding XR data of Funamori and Tsujii&. The ex- 
perimental 5*(g)'s show changes with increased pressure. 
The main peak grows in intensity and its position (q p ) 
increases monotonically, whereas the position of the sec- 
ond peak decreases between 8 and 14 GPa; the dis- 
tinctive shoulder at the high-q side of the main peak, 
shrinks smoothly and practically vanishes at 23 GPa. 
These changes are also reflected in g(r). The position 
of the main peak, (r p ), identified with the average near- 
est neighbor distance, decreases with pressure except for 
an increase between 8 and 14 GPa; the position of the 
second peak decreases monotonically with pressure. 

These features are displayed qualitatively in the calcu- 
lated S(qYs and g(r)'s although there are some quanti- 
tative discrepancies with the experimental data. Figure 
n shows that the OF-AIMD S(q) , s overestimate the in- 
tensity of the main peak and slightly underestimate the 
shoulder. Otherwise the positions of the peaks as well 
as the amplitudes of the subsequent oscillations are ac- 
counted for fairly well. A more detailed comparison with 
experiment is provided in Table [ri] which summarizes 
most of this structural information. This agreement with 
experiment is similar to that achieved in earlier orbital- 
free simulations 9 - and in KS-AIMDS*L£ calculations per- 
formed for for 1-Si near TP. 

Based on their experimental data for S(q) and g{r), 
Funamori and Tsujii& have argued that 1-Si undergoes 
a high pressure structural transformation between 8 and 
14 GPa. Whereas 1-Si contracts with pressure up to at 
least 8 GPa by reducing the bond length, as measured 
by r p , the increase in r p between 8 and 14 GPa suggests 
a structural change with an increase in the coordination 
number (CN). An estimate of the CN may be obtained 
by integrating the radial distribution function (RDF), 
4tit pi<?(r), up to the position of the first minimum, r m , 
in the RDF^iiS. The results from the calculated RDF 
in Table [ri] show that CN grows with compression but 
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FIG. 1: Static structure factor of 1-Si at different high pres- 
sures. Full circles: experimental X-ray diffraction data 16 . 
Continuous line: OF-AIMD simulations. 



with an abrupt increase from 8 to 14 GPa. Funamori 
and Tsujiiii obtained CN values also given in Table ITTIbv 
integrating their experimental RDF up to 3.1 A for all 
the states, and a similar growth between 8 and 14 GPa is 
seen. Had we used the same r m = 3.lA as Funamori and 
Tsuji did for calculating the CN, then the agreement with 
the "experimental" values would have been even better 
but, taking r m as the position of the first minimum is 
thought to be more soundly based. 

Values for the isothermal compresibility, kt , have been 
obtained from S(0) = piksTuT by using a least squares 
fit to calculated the S(q) for q- values up to 0.8 A -1 and 
extrapolating to q — > 0. Results are given in Tabic ITTT1 
Although no experimental results are available, the OF- 
AIMD calculation for 1-Si near the TP yielded kt = 3.0, 
which is rather close to the experimental valued of 2.8 
(in 10~ n m 2 Nw" 1 units). 

Further structural information is provided by higher 
order correlation functions such as the bond-angle distri- 
bution function, g3(9,r m ), where 9 is the angle between 
two vectors joining a reference particle with two neigh- 
boring particles at a distance less than r m . In a sim- 
ple liquid metal such as Al, the ga(0, r m ), has peaks at 
around 9 w 60° and 120°, which are close to those ex- 
pected for a local icosahedral arrangements 2 ^ (9 m 63.5° 
and 116.5°). In contrast, for 1-Si near the TP both the 
OF-AIMD2 and KS-AIMD&L& calculations for g 3 (0, r m ) 
have yielded two maxima centered around 9 « 60° and 
89°. This double-peak feature has been interpreted as a 
manifestation of tetrahedral bonding and higher coordi- 
nated atoms both contributing to the first coordination 
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FIG. 2: Bond-angle distribution function, <?3(0, r m ), of 1-Si at 
different high pressures (full line). The dotted and dashed 
lines stand respectively for 1-Si and 1-A1 at their respective 
TP. 



shell. In illustration, figure [21 shows the OF-AIMD re- 
sults for the g a (6, r m ) of 1-Si and 1-A1 near their TP'sASS. 
OF-AIMD results for compressed 1-Si are also included in 
the figure, and a gradual evolution with pressure towards 
the simple liquid metal distribution is seen. There is lit- 
tle change up to 8 GPa although the wide maximum at 
9 ss 89° has moved to slightly smaller 9 values which 
may indicate less tetrahedral bonding. But, as pressure 
increases from 8 to 14 GPa, there is a qualitative change 
in the g3(9,r m ) whose shape moves closer to that of the 
simple liquid metals, and at 23 GPa the positions of the 
maxima for 1-Si and 1-A1 are rather similar. 



TABLE II: Calculated values of q p (A -1 ), r v (A), r m (A) 
and coordination number (CN), for the different states. The 
numbers in parenthesis are the corresponding experimental 
data from Ref^- 



P(GPa) 


q P 


r p 


r m CN CN ° 


4 


2.61 (2.67) 


2.52 (2.46) 


3.03 6.6 (6.8) 6.9 


8 


2.65 (2.72) 


2.50 (2.42) 


3.07 7.2 (7.1) 7.3 


14 


2.78 (2.82) 


2.51 (2.46) 


3.20 9.6 (8.5) 8.9 


23 


2.84 (2.88) 


2.47 (2.43) 


3.28 11.0 (9.2) 9.5 



"Calculated by integrating the RDF up to r m = 3.1 A 
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FIG. 3: Dynamic structure factor S(q, ui) at several g-values 
(in A" 1 units), for 1-Si at 8 GPa. 



B. Dynamic properties. 

1. Collective dynamics. 

The intermediate scattering function, F(q,t), contains 
both spatial and temporal information on the collective 
dynamics of density fluctuations. It is defined as 




-iqR m (t+t ) 



m—1 




(1) 

Its frequency spectrum is the dynamic structure fac- 
tor, S(q,u), which has experimental relevance due to its 
connection with the scattered intensity in inelastic X- 
ray or neutron experiments. The calculated F(q, t) for 
compressed 1-Si exhibit an oscillatory behaviour up to 
q « (3/5) q p , with the amplitude diminishing for larger q- 
values. This oscillatory behaviour is typical of simple liq- 
uid metals found by either computer simulation22i2£i2£i2£ 
or from theoretical models 28 and gives rise to a well de- 
fined inelastic peak in S(q,uj). 

The S(q,uj), obtained by a time FT of F(q,t), exhibit 
for all the states, well defined sidepeaks which are indica- 
tive of collective density excitations. This is illustrated 
in figure |2| which shows calculated S(q,oj) for 1-Si at 8 
GPa for several q values. The general shape of S(q,u>) 
is qualitatively similar at equivalent q/q p values for all 
the compressed states, and no specific feature of S(q, lo) 
has been identified whose variation would mark the struc- 
tural transformation occurring somewhere between 8 and 
14 GPa. 

For all the states the sidepeaks in S(q,Lu) persists up 
to q w (3/5)q p , which is a feature shared by both l-SS 
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FIG. 4: Dispersion relation for the peak positions, uj m (q), 
from the calculated S(q, ui), for 1-Si at different high pressures. 
The figure also includes the calculated (open circles) (Ref£) 
and experimental (asterisks) (Ref^-) results for 1-Si near the 
triple point (T=1740 K). Dashed line: linear dispersion with 
the hydrodynamic sound velocity, v=3977 m/s., at the triple 
point. 



and the simple liquid metals near their TP'siSfli 2 ^ The 
dispersion relations, u m (q), of the density fluctuations 
have been obtained from the positions of these sidepeaks. 
They are plotted in figure El for the 8, 14 and 23 GPa 
states together with the calculated OF-AIMI32, results 
and experimental data— for 1-Si near its TP. The curves 
look qualitatively similar but there is a marked difference 
between uj m (q) for 8 GPa and less, and those for 14 and 
23 GPa. In addition, the slope of the dispersion gives a 
q-dependent adiabatic sound velocity, c s (q), which in the 
limit q — > reduces to the bulk adiabatic sound velocity, 
c s . This has been estimated by fitting a straightline to 
the low-g region of the u> m (q) 's and the results are given 
in Table ITTT1 The c s increase with pressure but a steeper 
rise from 8 to 14 GPa will be seen. 

The transverse current correlation function, J t (q,t), 
provides information on shear modes and is not directly 
related to any measurable quantity. It can only be ob- 
tained from either theoretical models or computer sim- 
ulations, but it is known that its shape evolves from a 
gaussian, in both q and t, at the free particle (q — > oo) 
limit, towards a gaussian in q and an exponential in t in 
the hydrodynamic limit (q — ► 0), i.e. 



Jt(q^0,t) = -L e -1 2 lW/™P> 



(2) 



where r\ is the shear viscosity coefficient, (3 = {ksT) 1 
and m is the atomic mass. In both small and large q 
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limits Jt(q,t) is always positive, but for intermediate q- 
values there is a more complicated behavior with well- 
defined oscillationsSS'SiiSS. Calculated Jt (q, t) for several 
q-values are shown in figures 15161 for 1-Si at 8 and 14 
GPa respectively. The most noteworthy effect on Jt(q, t) 
of an increasing pressure is reflected on the oscillations 
which have a smaller amplitude and last for appreciably 
shorter times at the lower presures. Its consequences are 
apparent in the frequency spectra, Jt(q,uj), which are 
plotted in the lower panels of figures 15161 For both 14 
and 23 GPa the Jt (<?, u>) exhibits an inelastic peak which 
appears at low q- values (« 0.45 A -1 ) and persists up to 
about q = 2.50 A -1 . However for 8 GPa the inelastic 
peaks appear for a appreciably smaller range (0.85 A" 1 
< q < 1.50 A -1 ) while for 4 GPa there are no inelastic 
peaks. This absence of peaks in J t (q,Lu) is also a feature 
of 1-Si near the TP, but is at variance with the behaviour 
of a large number of different liquids such as hard sphere 



systems**, Lennard- Jones liquids 



and simple liquid 



metals2SiSiiSl near melting, for which J t (q,t) oscillates 
and the associated Jt(q,uj), has an inelastic peak over 
some range of q- values. The inelastic peak in J t (q,u>) is 
associated with propagating shear waves which seem to 
be absent in 1-Si up to somewhere between 4 and 8 GPa. 
However, it must be noted that whereas in simple liquid 
metals near melting the shear waves last up to q « 3q p , 
in the case of 1-Si at 14 and 23 GPa they appear up to q 
« 0.9q p . 

The shear viscosity coefficient, 77, can be cal- 
culated from J t (q,t) using the memory function 
representation 20 ! 30 ! 31 



J~t(q,z) 



1 

[3m 



z -\ fj(q, z) 



(3) 



where the tilde denotes the Laplace transform, and 
fj(q, z) is a generalized shear viscosity coefficient. The 
J °° dtJt(q,t) when normalized gives (3m Jt{q,z = 0), 
from which fj(q, z = 0) are obtained and extrapolated to 
q = to give the usual shear viscosity coefficient, 77. Re- 
sults are given in Table llVl Although no comparison can 
be made with experimental results, the calculated values 
are considered reliable because the application of this ap- 
proach to 1-Si near its TP gave 77=0.75 ± 0.15 GPa • ps, 
in reasonable agreement with the available experimental 
data 3 ^, 77 exp =0.58-0.78 GPa • ps. It is noteworthy that 
once more r\ undergoes an abrupt change as the pressure 
increases from 8 to 14 GPa. 
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FIG. 5: (a) Transverse current correlation function Jt (q,t) at 
several g-values (in A -1 units), for 1-Si at 8 GPa. (b) Same 
for J t (q,uj). 



TABLE III: Calculated values of S(q — > 0), isothermal com- 
pressibility kt (in 10 -11 m 2 Nw" 1 units) and adiabatic sound 
velocity (c s ) for the different states. 



P(GPa) 


S(q 0) 


kt ( 10" 11 m 2 Nw' 1 ) 


c a (m/s) 


4 


0.0180 ± 0.003 


1.50 ± 0.3 


5100 


8 


0.0135 ± 0.003 


1.30 ± 0.3 


5400 


14 


0.0095 ± 0.003 


0.94 ± 0.3 


6300 


23 


0.0085 ± 0.003 


0.68 ± 0.3 


6750 
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FIG. 6: (a) Transverse current correlation function Jt (q,t) at 
several (/-values (in A -1 units), for 1-Si at 14 GPa. (b) Same 
fo J t (q,w). 



2. Single-particle dynamics. 

Information about the single-particle properties is con- 
tained in the self-intermediate scattering function 



F s (q,t) = _(^ e -^i(*+*o) e i ^( t0 )) (4) 

3=1 

and its frequency spectrum, the self-dynamic structure 
factor, S s (q,uj), which is related to the incoherent part 
of the total intensity scattered in an INS experiment. 
The OF-AIMD results for the F s (q, t) presented in figure 
0for two states display the usual monotonic decay with 
time, and comparison of different states shows that at 
similar q/q p values F s (q,t) decays slower with increasing 
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FIG. 7: Self intermediate scattering functions, F s (q,t), at 
several g-values (in A -1 units), for 1-Si at (a) 8 GPa and 
(b) 14 GPa. The key for the symbols is the same as in the 
previous figure. 



pressure, with the 14 and 23 GPa states behaving very 
much like the liquid simple metals22i2i2& near their triple 
points. The different rates of decay can be related to the 
differences in the self-diffusion coefficients. 

Closely related to the F s (g, t) is the velocity autocorre- 
lation function (VACF) of a tagged ion in the fluid, Z(t), 
which can be obtained as the q — > limit of the first-order 
memory function of the F s (q,t), but more conveniently, 
from its definition 



z(t) = mwmM) (5) 

Figures T8I9I show results for Z(t) and for its power spec- 
trum Z(u>). The overall shape of Z(t) changes little from 
the TP up to 8 GPa closely, but as the pressure is fur- 
ther increased changes occur in the range and amplitude 
of the oscillations until at 23 GPa the shape is very simi- 
lar to that of the simple liquid metals at their TP. These 
changes can be explained in terms of the so-called " cage" 
effect due to backscattering from the shell of nearest 
neighbors reversing the initial velocity of a tagged ion 
and driving a deeper first minimum. This is consistent 
with the results for the static structure summarized in 
table ITU which shows an open structure up to 8 GPa but 
a marked increase in the coordination number at higher 
pressures. 

The power spectra which are plotted in figure [5] also 
show significant changes between 8 and 14 GPa. Z(u>) 
evolves with pressure from a shape resembling 1-Si near 
the TP towards a liquid simple metal shape with a 
low frequency peak and a higher frequency peak (or 
shoulder The shoulder at ui « 40 ps -1 , present 
at all the pressures, has been related to vibrational rem- 
nants in the liquid of the covalent bonding 6 . Below 8 
GPa, low frequency diffusive modes are present. 

The self-diffusion coefficient, D, is readily obtained 
from either the time integral of Z(t) or from the slope 
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FIG. 8: Full line: Normalized velocity autocorrelation func- 
tion Z(t) for 1-Si at different pressures. The dashed and dash- 
dotted lines stand for the respective longitudinal, Zi(t), and 
transverse, Z t (t), components as defined in equation 0. The 
dotted line depicts the result for 1-Si at triple point (T =1740 
K). 

of the mean square displacement 5R 2 (t) = (\Ri(t) — 
Ri(0)\ 2 ) of a tagged ion in the fluid, as follows 

1 f°° 

D = — / Z(t)dt ; D = lim SR 2 (t)/6t (6) 
pm J t—>oo 

Both routes for D lead to practically the same value, 
and the results are given in Table CSI The decreasing 
values of D with increasing pressure is due to the grow- 
ing importance of backscattering. No experimental re- 
sults are available for the diffusion coefficients of 1-Si, but 
confidence in the results may be taken from the agree- 
ment between the OF-AIMD result for 1-Si near the TP: 
-Dof— aimd = 2.28 A 2 /ps, and the estimates from KS- 
AIMD calculations of Stich et afr£: -Dks-aimd = 2.02 
A 2 /ps£ which slightly increased to 2.4 A 2 /psi£ when the 
number of particles was augmented to 350 particles. An- 
other KS-AIMD study by Chelikowsky et aim has yielded 
-Dks-aimd =1-90 A 2 /ps. The results for D at 4 and 
8 GPa are similar to the value at the 1-Si triple point 
whereas the results for 14 and 23 GPa are closer to those 
for the liquid simple metals near their triple points22i2&2i. 
This change in D with pressure explains the different de- 
cay rates found in the F s (q, t) which decayed much faster 
at the lower pressures. Recalling the accurate gaussian 
approximatio n 24 ! 29 , F s (q,t) — exp[~q 2 SR 2 (t) /6], it will 
be seen that a greater D implies a greater 5R 2 (t) and 
therefore a faster decay of F s (q,t). 



FIG. 9: Power spectrum, Z(uf), for 1-Si at different pressures. 
The dashed and dash-dotted lines stand for the respective 
longitudinal, Zi(u>), and transverse, Zt(oj), components. The 
dotted line depicts the result for 1-Si at triple point (T =1740 
K). 



TABLE IV: Calculated values of the self-diffusion (D) and 
shear viscosity rj (in GPa ps) for the different states. 



P(GPa) 


D (A 2 /ps) 


77 (GPa ps) 


4 


1.82 ± 0.05 


0.77 ± 0.10 


8 


1.33 ± 0.05 


0.84 ± 0.10 


14 


0.70 ± 0.03 


1.47 ± 0.15 


23 


0.70 ± 0.03 


1.55 ± 0.15 



The self-diffusion coefficient, D, of a macroscopic parti- 
cle of diameter d undergoing Brownian motion in a liquid 
of viscosity 77 is related to 77 through the Stokes-Einstein 
(SE) relation 77 D = k^T '/2nd. This relation has often 
been used on an atomic scale to estimate 77 by identifying 
d with the position, r p , of the main peak in g(r). Using 
the D values for 4, 8 ,14 and 23 GPa the relation yields 
77=0.72, 0.82, 1.36 and 1.59 GPa-ps respectively, values 
rather close to the earlier OF-AIMD estimates. 

Gaskell and Miller—, have used mode-coupling the- 
ory to develop a representation of the normalized VACF 
which has been used to interpret MD data in various 
fluido 35 ! 36 ! 37 , and which sheds light on 1-Si. Within this 
approach 

Z(t) «-^3 f dqf(q)[J l (q,t) + 2J t (q,t)]F s (q,t) 
= Zi(t) + Z t (t) 

(7) 

where Ji(q, t) and Jt(q, t) are the normalized longitudinal 
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and transverse current correlation functions and f(q) is 



2.0 



/(g) 



3 ii(ag) 

Pi aq 



(8) 



with the spherical Bessel function of order one, pi is 
the ion number density and a = (3/47rpi) 1 / 3 is the radius 
per ion. Substitution into equation J2J of the OF-AIMD 
results for J;(g, t), J t (q,t) and F s (q,t) allows identifica- 
tion of longitudinal and transverse current contributions: 
Zi(t) and Z t (t) respectively. The two contributions are 
plotted in figure |H1 which shows that the oscillatory be- 
haviour in the Z(t) is due to Zi(t), but again the step 
from 8 to 14 GPa changes the shape of both contribu- 
tions. Up to 8 GPa, Z t (t) remains positive for all times 
as a result of the positive nature of Jt(q,t) (see figures 
151611 and determines the long time behaviour of the Z(t), 
however, from 14 GPa on, Z t (t) develops a shallow and 
broad negative minimum centered at rather long times ss 
0.15 ps. On the other hand, Z[(i) accounts for most of the 
backscattcring effect. With higher pressure the first min- 
imum sharpens and moves to shorter times and the oscil- 
lations extend further, which are results of the increasing 
role of the "cage" effect. At 14 and 23 GPa, both com- 
ponents are similar in shape to their liquid simple metals 
counterparts^ with both oscillating about zero and with 
Zi(t) controlling the large t behaviour of Z(t). Finally, 
notice that the development of the deep minimum in the 
Z(t) is mainly due to the rapid decay of Z t (t) with in- 
creasing pressure. 

The longitudinal and transverse components of the 
power spectrum, Z(oS), are shown in figure^l The spec- 
trum at small to is dominated by Z t {uj), and, conse- 
quently, the diffusion constant D oc Z(u> — 0) is com- 
pletely determined by the transverse component. For 4 
and 8 GPa, Z t {u>) decreases monotonously but at 14 and 
23 GPa the value at zero frequency has dropped and a 
low- frequency peak has developed. Note that Z t (ui) has 
no maximum for 4 and 8 GPa which are the states where 
the Jt(q,uj) shows, either no inelastic peaks (4 GPa) or 
they exist for a small range (8 GPa). The longitudinal 
component Zi(lu) always exhibits a peak whose position 
increases slightly with increasing pressure. This peak is 
responsible for the shoulder in the total Z(u>) for the 4 
and 8 GPa states, as well as for the high-frequency peak 
for the 14 and 23 GPa states. 

By a time FT of the F s (q,t) we obtain its frequency 
spectrum, S s (q,uj), which is known as the self-dynamic 
structure factor and is related to the incoherent part of 
the measured INS cross-section. For all q- values S s (q, ui) 
decays monotonically as a function of frequency and can 
be characterized in terms of the peak value, S s (q, lo = 0) 
and the HWHM, LUi/ 2 (q)- These parameters are fre- 
quently reported normalized with respect to the hydro- 
dynamic (q — > ) limit, by introducing the dimension- 
less quantities S(q) = irq 2 DS s (q, u> — 0) and A(q) — 
uji/2(q)/q 2 D, where u)i/2(q) / q 2 can be interpreted as an 
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q(A ) 



FIG. 10: Normalized HWHM of S s (q, w), relative to its value 
at the hydrodynamic limit, for 1-Si at 8 and 14 GPa. Con- 
tinuous line: OF-AIMD results. Asterisks: Mode-coupling 
theory. Dashed line: free-particle limit. 
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FIG. 11: Same as before, but for the normalized peak value 
S s (q,uJ = 0), relative to its value at the hydrodynamic limit. 



effective g-dependent diffusion coefficient D(q). For a 
simple liquid near its TP, A(q) usually oscillates whereas 
in a dense gas it decreases monotonically from unity at q 
= to the 1/q behaviour at large g 24 i 25 i 29 . Fieures HOIllI 
depict the OF-AIMD results for A(q) and S(g) for 8 and 
14 GPa, as this is the pressure range where both magni- 
tudes undergo a substantive change. The obtained A(q) 
for 14 GPa shows an oscillatory shape with a minimum 
located at q w q p which can be traced back to structural 
features ("cage" effect) which somewhat hinder the mo- 
tion of the ions and becomes more effective at q w q p 
where the wavelength is comparable to the size of the 
cage. Conversely, that for 8 GPa resembles the situation 
of a dense gas where the " cage" effect becomes negligible 
leading to a net reduction of the diffusion coefficient! 2 ^^ 
and its associated £(<?), stands very close to that of the 
dense gas. 

An additional check on the reliability of these results 
may be provided by the MC theory22i2P. which has already 
shown its capability to describe to experimental data for 
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A(q) and £(<?) in simple liquid metals^Siii at q < q p . 
Specifically, the MC theory avers that at low-g values 

A(q) = 1 + H (S)q/q* (9) 

where q* = 16nmpi(3D 2 , 5 = D/(D + r\jmpi) and H(S) 
and G(5~ 1 ) are given in Reference^. 

The first term in equations 10 stands for the hydro- 
dynamic result whereas the second one accounts for the 
coupling of mass diffusion and collective modes. Calcu- 
lated values of D and r\ have been used with equations 
El to obtain the points in figures 11011 II For q < q p we 
observe that the MC theory fairly accounts for the OF- 
AIMD results, with an accuracy comparable to what has 
already been achieved in other liquid metals^. Conse- 
quently, the present results show the ability of the MC 
theory to describe the single particle dynamics (and pre- 
sumably the collective dynamics too) in liquid systems 
encompassing a range of bonding and structure as that 
displayed by the compressed 1-Si . 

IV. CONCLUSIONS. 

Several static and dynamic properties of 1-Si at four 
high-pressure thermodynamic states have been investi- 
gated using orbital free ab initio molecular dynamics 
combined with a first-principles local pseudopotential. 

The study was motivated by experimental findings^ of 
significant structural changes in 1-Si when the pressure is 
increased from 4 to 23 GPa. The obtained results for the 
static structure qualitatively follow those trends unveiled 
by the experiment, namely the increase of the intensity 
and the position of the S(q)'s main peak, along with a 
progressive vanishing of its shoulder. Other parameters 
such as the coordination number, isothermal compress- 
ibility, and the shape of the bond-angle distibution func- 
tion provide further insight into the changes. Overall, 
apart from a contraction with increasing pressure, the 
static structures of 1-Si at the TP, and at 4 and 8 Gpa 
are very similar. Above 8 GPa the system transforms to 
a denser more close packed structure typical of a liquid 



simple metal, with most change taking place between 8 
and 14 GPa. 

The structural changes are also reflected in several dy- 
namical properties. The calculated dynamic structure 
factors, S(q,u)), show collective density excitations over 
similar wavelength ranges, namely up to q « (3/5)g p , as 
those found for simple liquid metals at their TP. These 
density excitations are sound waves whose velocity in- 
creases with pressure, most steeply between 8 and 14 
GPa. The dispersion relations of the excitations divide 
into two groups, one for 1-Si at its TP, 4 GPa and 8 GPa 
and another group for 1-Si at 14 and 23 GPa. 

The transverse current correlation also show evidence 
of the structural changes. Below 4 GPa, its frequency 
spectra lack inelastic peaks indicating the absence of 
shear waves, but at 8 GPa clear inelastic peaks are al- 
ready evident. 

The calculated self-diffusion and shear viscosity trans- 
port coefficients are also affected by the structural 
changes occurring between 8 and 14 GPa. These trans- 
port coefficients cannot be compared with experiment, 
but confidence in the calculated values is given by the 
good agreement with experimental values and/or other 
ab-initio results for 1-Si and its TP. 

Finally, we remark that the present results for the 
static and dynamic properties of compressed 1-Si under- 
score the capability of the OF-AIMD method to tackle 
liquid systems encompassing a range of bonding and 
structure which evolves from mild remnants of covalent 
bonding to a metallic one. Moreover, further improve- 
ments in the present ab initio method are still possi- 
ble and they necessarily will be focused on developing 
more accurate electron kinetic energy functionals and lo- 
cal ionic pseudopotentials. 
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